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The phase- field-crystal model for liquid crystals is solved numerically in two spatial dimensions. 
This model is formulated with three position-dependent order parameters, namely the reduced trans- 
lational density, the local nematic order parameter, and the mean local direction of the orientations. 
The equilibrium free-energy functional involves local powers of the order parameters up to fourth 
order, gradients of the order parameters up to fourth order, and different couplings between the 
order parameters. The stable phases of the equilibrium free-energy functional are calculated for 
various coupling parameters. Among the stable liquid-crystalline states are the isotropic, nematic, 
columnar, smectic A, and plastic crystalline phases. The plastic crystals can have triangular, square, 
and honeycomb lattices and exhibit orientational patterns with a complex topology involving a sub- 
lattice with topological defects. Phase diagrams were obtained by numerical minimization of the 
free-energy functional. Their main features are qualitatively in line with much simpler one-mode 
approximations for the order parameters. 
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I. INTRODUCTION 

Crystallization and melting processes (see e.g. [TJ [2]) 
can be efficiently modelled within the so-called phase- 
field-crystal (PFC) model [3HS] which is basically a 
Landau-type theory with a conserved position-dependent 
order parameter which is gradient-expanded up to fourth 
order. Stable periodic oscillations in the order are in- 
terpreted as crystalline density fields. The PFC model 
was applied to various situations including fluid-crystal 
interfaces [6j [7], crystal growth into a supercooled liq- 
uid [8], grain boundaries [9j [10], and the Asaro-Tiller- 
Grinfeld instability [TTHT4] . Later on, the PFC model 
was systematically derived from density functional the- 
ory [UEJ[8l[T5j[T6] which provides a microscopic approach 
to freezing and melting in equilibrium [T7H20] using a gra- 
dient expansion in terms of density modulations p~6j [2T1 
[23] . Moreover, dynamical density functional theory for 
the non-equilibrium dynamics of Brownian systems [24] - 
[26] can be used to derive the dynamics of the PFC model 
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Recently, the PFC model was generalized to 
anisotropic particles by using the appropriate generaliza- 
tions of density functional theory to orientational degrees 
of freedom [20], [27j [28] and to orientational Brownian dy- 
namics [29 , 30 . This was done both in two 31 and three 
[32] spatial dimensions [33] . The resulting PFC models 
are valid for liquid crystals composed of apolar particles 
and describe also liquid-crystalline phases with a uniax- 
ial orientation distribution. The theory is formulated in 
terms of three order parameter fields, namely the reduced 
translational density, the local nematic order parameter, 
and the orientational direction. It includes gradients up 
to fourth order in the reduced translational density field 
and up to second order in the remaining orientational 
order parameters. While the original PFC model has 



two free parameters, the liquid-crystalline PFC model 
in two dimensions [31] has five independent couplings. 
This widely opens the parameter space for the occur- 
rence of several liquid-crystalline phases including ne- 
matic, columnar, smectic A, plastic crystalline, and ori- 
entationally ordered crystalline phases [34, 35J. However, 
no numerical calculation of the PFC model has been pre- 
sented for the stability of these different liquid-crystalline 
phases, yet. 

In this paper, we start from the PFC model for liq- 
uid crystals in two spatial dimensions proposed in ref- 
erence [31] and determine the stable liquid-crystalline 
phases numerically for special coupling-parameter com- 
binations. Among the stable liquid-crystalline phases are 
the isotropic, nematic, stripe, columnar, smectic A, and 
plastic triangular crystalline phase. The latter exhibit 
complex orientational patterns with a sublattice of topo- 
logical defects. For stronger translational-orientational 
couplings, a plastic crystal with fourfold square symme- 
try is getting stable. This does not occur for the tradi- 
tional PFC model [3] [4] [36] but can be induced for other 
Ginzburg-Landau functionals [37] . Also a plastic honey- 
comb crystal was observed. We found second-order phase 
transitions for the isotropic-nematic phase transition and 
phase transitions of first order for all other phase trans- 
formations. 

The paper is organized as follows: S ec. [TT] is addressed 
to the PFC model for liquid crystals. In |II A[ we recall the 
free-energy functional for two-dimensional liquid crystals 
and introduce a simpler form of this functional, which is 
obtained by an alternative choice of the order-parameter 
fields. Two different methods for the minimization of this 



functional are presented in II B Then, in Sec. |III[ results 
for these two methods are shown. We finally conclude in 
Sec. El 
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II. PHASE-FIELD-CRYSTAL MODEL 

A. Free-energy functional 

We start our investigation from the PFC model for 
liquid-crystalline phases, that was recently proposed in 
reference [3T. It can be rescaled to the following form: 
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This functional approximates the free energy of a set of 
apolar uniaxial particles with the three space-dependent 
order-parameter fields i/ji(x), ^O?), and &o0*0- The field 
(5) is the reduced orientationally averaged one-particle 
density obtained from the general one-particle density 
p(x,u) which corresponds to the probability density to 
find a particle with the orientation u at the position x as 
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du (p(x, u) — p) . 



Here, S± is the unit circle, p is the mean density 
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and the unit vector u denotes a certain orientation in 
the two-dimensional space. In comparison to the gen- 
eral one-particle density p(x,u), the reduced and ori- 
entationally averaged density ipi(x) only describes the 
space-dependent deviations of the one-particle density 
from the reference density p. Its spatial average van- 
ishes by construction. The orientational dependence of 
the free-energy is taken into account by the nematic order 
parameter ^(x) 
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and the nematic director uo(x), which is the eigenvec- 
tor associated with the largest eigenvalue of the second- 
order traceless nematic tensor Q(x) [38]. The nematic 
director uq(x) is a unit vector and can be parametrized 
by the orientation field (fo(x) in two spatial dimensions: 
uq(x) = (cos((fo(x)) ^ sm(cpo(x))). Also the amount of lo- 
cal ordering, ^2 0*0? contributes to the free energy. In 
the integrand of functional 0, the first four terms ap- 
proximate the free-energy density of an ideal rotator gas. 
Also the terms ip 2 (x) and ^ 2 {x) appear in the ideal ro- 
tator entropy, but since there are corresponding terms 



in the excess free energy for anisotropic particles, the 
different contributions to these terms were combined to 
B\ipi(x) and Difj^x) in Eq. 0. Beside the mentioned 
polynomial terms in the order-parameter fields, also their 
gradients contribute to the free energy. The amount of 
their contribution is controlled by the parameters £? x , E, 
and F in the free-energy functional. Contributions of the 
gradients and curvature of the translational density field 
ipi(x) appear in the term proportional to £? x . The cur- 
vatures of the space-dependent nematic order ^2 (x) and 
of the nematic director field iio(x) are taken into account 
in the term proportional to the parameter E. The last 
term in the free-energy functional is scaled by the param- 
eter F and contains the couplings between the gradients 
of ipi(x) and ^(x) as well as a sum, which couples the 
Hessian of ipi(x) with the director field uo(x) and with 
the nematic order ^(x). In Eq. Q, the components 
of the nematic director are denoted as ui = (u(x))i, for 
abbreviation. 

It is possible to transform the functional (TJ to a physi- 
cally less intuitive but simpler equation by defining a new 
complex order-parameter U(x) = ^2(x)e l2ipo ^ x \ The cor- 
responding free-energy functional 
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with U*(x) denoting the complex conjugate of U(x) are 
helpful for more efficient numerical simulations. This al- 
ternative choice for the order-parameter fields is equiv- 
alent to the usage of a not normalized vector field in- 
stead of iio(x) and the modulus ^2 0*0? as it is done in 
Ginzburg-Landau theory [38] . 

Before we started to minimize the free-energy func- 
tional numerically, we could already derive some static 
properties of this functional. Since the parameter D 
controls the contribution of the nematic order parameter 
^2 (x) , we expect the nematic phase to be stable for large 
negative values of D. In the opposite case, if D is large 
enough and positive, the term 0*0 +^2 0*0 /256 domi- 
nates the free energy and only phases with ^(x) = can 
be stable. Crystalline phases with a non- vanishing ne- 
matic order might only appear in a region around D = 0, 
therefore. We also know that the difference B\ — B x has 
a big influence on the translational density field ^i(x). 
If the parameter B\ is large and positive, variations of 
the translational density field enlarge the free energy. 
Similarly, gradients of the translational density field en- 
large the free energy for large and negative values of B x . 
Therefore, phases without any density modulations, i.e. 
the isotropic and the nematic phase, are preferred for pos- 
itive values of the difference B\ — B x . All other phases 



3 



isotropic 



nematic 




plastic triangular crystal 1 



oo 



•I 



r _ i 



plastic honeycomb crystal 



v 



X 10 





10 15 5 

columnar/smectic A 





10 12 2 4 6 x 8 10 12 

plastic triangular crystal 2 




^2 (X) 



FIG. 1. Stable liquid crystalline phases. The contour plots show the order-parameter fields i/ji(x) and 1^2 (x) for the isotropic 
and nematic phase, the stripe phase and columnar/smectic A phase, two plastic triangular crystals with different orientational 
ordering, and a plastic honeycomb crystal as well as a plastic square crystal. The black lines in the plots of the second and 
fourth column represent the director field uo(x). In the cases with ip2(x) = it is not shown, because it is not defined. The 
parameters are = 3.5, E — 1, and F — for the stripe phase and the plastic triangular crystal 1 and — 3.5, E — 0.1, 
and F = 1 for all other phases. 



with a periodic translational density field are preferred 
for negative values of this difference. Furthermore, there 
is a symmetry concerning the reversal of the sign of the 
parameter F in the free-energy functional. From Eq. 
we know that the free-energy functional is invariant un- 
der a simultaneous change of the signs of the parameter 
F and the nematic order-parameter field ^{x). Due to 
this symmetry, we can assume F ^ in the following. 



B. Numerical minimization 

1. Free numerical minimization 

In order to find the stable phases in the PFC model, 
we minimized the free-energy functional ([T]) using the 
steepest descend method for fixed parameters B\, £? x , D, 
E, and F. In this method, which is based on the pseudo- 
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FIG. 2. Topological defects in three different plastic liquid 
crystals (schematic). The defects coincide with the maxima 
(red discs) and minima (cyan discs) of the translational den- 
sity field ipi (x) . The symbols in the plots represent the follow- 
ing defects: vortices (a) with the topological winding number 
m = 1, disclinations (b) with m — —1/2, sources/sinks (c) 
with m = 1, and hyperbolic points (d) with m = — 1. 



dynamical equations 

di/ji(x,t) ST^^x, t), U(x, t)} 

dt Sipi (x, i) 

dU(x, t) ST^jx^.U^t)] 
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the system evolves towards a local minimum. Here, A 
is a Lagrange multiplier that guarantees that the spatial 
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FIG. 3. Phase diagrams calculated by full numerical mini- 
mization for the parameters — 3.5 and E — 0.1. The 
relevant liquid crystalline phases are: isotropic (blue), ne- 
matic (green), stripes (yellow), columnar/smectic A (C/SA, 
orange), plastic triangular crystals (magenta), plastic honey- 
comb crystal (PHC, dark purple), and plastic square crystal 
(PSC, red). 



average of the field ^i(x) is zero. The two equations were 
discretized using a finite-difference scheme and solved on 
a grid with 32 x 32 cells. We varied the length of these 
cells to minimize the free-energy functional also with re- 
spect to the lattice constant of the periodic phases and 
used a set of different phases as initial conditions in order 
to find the global minimum. 



2. One-mode approximation 

As a semi-analytical approach the one-mode approxi- 
mation consists of periodic approximations for the order- 
parameter fields i/ji(x), "02 0?) ? an d (po(x) and reduces the 
PFC model to the lowest Fourier modes. We used the fol- 
lowing parametrizations with the minimization parame- 
ters A\, £?o, and k: 



isotropic phase: 



i/j 2 (x) 



0, 
0, 



nematic phase: 
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• stripe phase: 



i/ji(x) = A x cos(ky) , 

^2(5) = 0, 

• columnar phase: 

^i(f) = A x cos(ky) , 
ip2(x) = B + Bx cos(ky) , 
cp (x) , 

• smectic A phase: 

ipx(x) = Ax cos(ky) , 
foix) = B + Bx cos(ky) , 

<po(x) 
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triangular crystalline phase: 

i/>i(x) = Ax (cos (~Y kx ) cos i^v) 

cos(ky)\ 

^ 2 (f) = B + #i ^cos (-^fcz) cos (^2/) 
cos(ky) 



(12) 



<p (£) </>o , 
• square crystalline phase: 

ipi(x) = Ax ( cos(/cx) + cos(%)) , 
^2 0*0 — + #i cos(/cx) cos(ky) , 
<Po(#) • 

Here, x and ?/ are the first and the second Cartesian coor- 
dinate, respectively. The constant angle </>o can be set to 
zero, because the free energy does not depend on it. Fur- 
ther, due to equivalent free energies, we do not need to 
distinguish between the columnar phase and the smectic 
A phase and call them columnar /smectic A phase in the 
following. Note, that there is no additive offset term Ao 
for the density variations ipx (x) . The minimization of the 
free energy in the context of the one-mode approximation 
was performed for fixed parameters B\, £? x , D, E, and F. 
It is possible to minimize the free energy for the nematic 
phase and for all phases with a vanishing nematic order 
ip2(x) analytically. The more complicated free energies 
of the remaining phases were minimized numerically. 



III. NUMERICAL RESULTS 
A. Free numerical minimization 

Apart from the fully isotropic phase (^x(x) = 0, 
ip2(%) — 0), which appears for B\ > £? x and D > 0, sev- 
eral other phases were found to minimize the free-energy 



(see figure [T]). As expected, for negative and large D, 
a nematic phase was found. In the columnar/smectic A 
(9) phase, the system has positional ordering in one direction 
while it is isotropic perpendicular to this direction. The 
ip2(x) field has a similar profile to the ipx(x) field with 
maxima of these two fields at the same positions. The 
director field uq(x) is oriented parallel to the gradient 

(10) of the field i/ji(x), V^i^), in the regions with negative 
Vi^i(x), while for positive V^i^) it is perpendicular 
to the gradient. A similar flipping of the orient ational 
field from perpendicular to parallel to the stripe direc- 
tion was identified as transverse intralayer order in the 
three-dimensional smectic A phase of hard spherocylin- 

(11) ders [39]. 

Four plastic crystals with different symmetries were 
found. The first two phases are plastic triangular crys- 
tals with a vanishing and a non- vanishing nematic order 
parameter. Here, the first crystal with ^(x) = is a 
special degenerate case of the second. The third plastic 
crystalline structure involves a honeycomb lattice. As a 
fourth case, there is a plastic crystal with square symme- 
try. For all plastic crystals, ip2(%) vanishes both at the 
maxima and minima of ipi(x). 

The director fields uq(x) for the different plastic crys- 
talline phases exhibit quite different topologies. While 
the director field is not defined for the plastic triangu- 
lar crystal 1 with a vanishing field ^2 (a?)? it possesses in 
general topological defects at positions where the field 
ip2(x) vanishes. This guarantees that there is a finite 
core energy of the topological defect [40]. The topologi- 
cal defects form another lattice with more lattice points 
than that given by the maxima of the field i/)i(x) since 
(13) there are additional interstitial topological defects at the 
minima of ^\(x). The lattices of topological defects are 
schematically schown in figure [2] For the plastic triangu- 
lar crystal 2 and for the plastic honeycomb crystal, the 
associated defect crystal is triangular albeit with a lat- 
tice constant that is a factor of 1/a/3 smaller than the 
original one. Likewise, for the plastic square crystal, the 
defect lattice is a square lattice with a lattice constant 
reduced by a factor l/y/2. Topological defects in liq- 
uid crystals can be classified according to the winding 
number of their director field [38 , 41 . In our case, three 
types of point defects occur: vortices with the topological 
winding number m = 1, sources/sinks with m = 1, and 
hyperbolic points with m = — 1. Furthermore, we found 
disclination line defects with m = —1/2. Vortices and 
disclination lines occur in the plastic triangular crystal 
2, which is schematically drawn in the first plot in figure 
|2j In the plastic honeycomb crystal, disclination lines 
arise together with sources/sinks (second plot in figure 
|2j, while vortices and hyperbolic points are found in the 
plastic square crystal (last plot in figure [2]). The sum of 
the topological winding numbers of all topological defects 
in a unit cell vanishes for all plastic crystals. 

The director fields of all crystalline phases that we 
found are periodic and vanish when they are averaged 
in space. Therefore, these crystals are identified as being 
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plastic. Orientationally ordered crystalline phases were 
not found in the parameter range we explored. 

Parameters for which the mentioned phases are stable 
follow from the phase diagrams in figure (3) We chose 
1.5 ^ B x < 4.5, B x = 3.5, -1 ^ D < 1, E = 0.1, and 
F = 0,0.01,0.1,1 for these phase diagrams. We chose 
B x constant and varied B\ so that in the absence of the 
orientational degrees of freedom we obtain all the phases 
of the original PFC model. The parameter D was varied 
from —1 because for this value the energy of the nematic 
phase was significantly smaller than the energy of any 
other phase of the original PFC model. The parameters 
E and F were selected from the regions where we found 
the richest phase diagrams. For F = we have a degen- 
erate case, where ip2(%) — for D > and E > 0. In 
the latter case we observed a stripe phase and the plastic 
triangular crystal 1. These phases are replaced by the 
columnar /smectic A phase and by the plastic triangular 
crystal 2 if F becomes positive. The richest phase dia- 
gram with six different phases was obtained for F = 1. 
The phase transition between the isotropic and the ne- 
matic phase turned out to be continuous, while all other 
phase transitions are discontinuous. This result agrees 
with the fact that the PFC model reduces to the Landau- 
de Gennes model for the isotropic-nematic phase transi- 
tion, which describes this phase transition as continuous. 
On the other hand, symmetries are broken for the re- 
maining phase transitions so that they are discontinuous. 



B. One-mode approximation 

For the global minimization of the free energy we used 
a random search routine in the four-dimensional param- 
eter space combined with a local minimization by a New- 
ton method. Figure [4] shows the order-parameter fields 
ipi (x) and ^2 (x) for different phases that are realized for 
particular combinations of the parameters of the PFC 
model in the one-mode approximation. To calculate 
phase diagrams that correspond to the phase diagrams 
presented in Subsec. |III A[ we fixed the parameters £? x , 
E, and F to the same values as in the previous sub- 
section, varied the remaining parameters of the rescaled 
free-energy functional over the intervals — 1 ^ D ^ 1 
and 1.5 ^ B\ ^ 4.5, and minimized the free-energy for 
each one-mode approach and for each parameter combi- 
nation. The resulting phase diagrams are shown in figure 
[5] A view at the behavior of the order parameters near 
to phase transitions confirmed the results of the free nu- 
merical minimization regarding the order of the phase 
transitions. 



method to determine the phase diagram for the PFC 
model. A speed-up factor of 100 can easily be achieved. 
This allows the fast calculation of phase diagrams for 
various parameter combinations as well as single phase 
diagrams with a rather high resolution, like those that 
are shown in figure [5) The drawback, however, is that 
the one-mode approximation also gives rise to deviations 
of the phase diagrams from the exact ones in figure [3J 
there is a lack of some crystalline phases for F = 1. The 
absence of these crystalline phases is based on an im- 
proper consideration of the orientation field ipo(x) in the 
one-mode approximation and the last term in the free- 
energy functional ([!]) which is relevant for non- vanishing 
values of the parameter F. Therefore, the plastic square 
crystal as well as the plastic honeycomb crystal cannot 
be observed in the one-mode approximation. In addition, 
the plastic triangular crystal 1 with ^(x) = appears 
for F > instead of the plastic triangular crystal 2. This 
is not surprising since we used a constant angle 0o for the 
orientation field in our ansatz for the one-mode approx- 
imation to parametrize plastic crystals without a global 
orientation. The constant angle </>o is not a good approx- 
imation for the actual orientation field with a vanishing 
mean value so that the minimum of the free energy is 
reached for ^(x) = 0. A further difference between the 
phase diagrams in figure [5] and those in figure [3] is an is- 
land of the columnar /smectic A phase near to the lower 
end of the phase transition line between the isotropic and 
the nematic phase for F = 0.1. 

To check that the mentioned disadvantages of the 
phase diagrams for the one-mode approximation really 
arise from the improper approximation of the orienta- 
tion field <po(%), we performed a Fourier analysis of the 
numerical results that are shown in figure [I] This Fourier 
analysis exhibited that the first Fourier mode is dom- 
inant. Only for the plastic square crystal and for the 
plastic honeycomb crystal there appears a contribution 
of the second mode that might be relevant. 

Despite the discussed disadvantages of the one-mode 
approximation, it proved to be an useful method in order 
to calculate phase diagrams for a given PFC model with 
low computational effort. The one- mode approximation 
is also useful to explore the phase diagram for suitable 
parameters or for a large number of parameter sets and to 
find interesting regions in a high-dimensional parameter 
space that are worth to be investigated more precise by 
a much more expensive direct minimization of the free- 
energy functional. 



IV. CONCLUSIONS 



C. Comparison and discussion 

In comparison to the full numerical minimization of 
the free-energy functional, the one-mode approximation 
appears to provide an approximative but much faster 



In conclusion, the phase-field crystal model for liquid 
crystals was solved numerically in two spatial dimensions 
and different stable phases were found in the parameter 
space of the model. These include isotropic, nematic, 
columnar, smectic A, and plastic crystalline phases. The 
latter can possess a periodic triangular, square or hon- 
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eycomb structure with a complex orientation field. The 
phase-diagrams were obtained by a numerical minimiza- 
tion of the free-energy functional but are - in most but 
not all cases - qualitatively in line with much simpler 
one- mode approximations for the order parameter fields. 
We hope that this work will inspire more simulations and 
experiments to observe the unexpected predicted phases 
(like the square and honeycomb crystals) in two dimen- 
sions [42H46] . One can think about the realization of 
molecular liquid crystals as well as concentrated solutions 
of anisotropic colloidal particles [47[ 08] or anisotropic 
mesoscopic dust particles in a plasma [49, 50 . In partic- 
ular, it would be challenging to explore the orientation 
field and the associated defect lattice in plastic crystals 
either by simulations [51 or by experiments [52) [53]. 

An important further step is to calculate the full liquid- 
crystalline phase diagram for a given interparticle poten- 
tial as a function of the real thermodynamic variables, 
namely temperature and number density. This phase di- 
agram is known from computer simulations and from the- 
ory for hard spherocylinders [35j [54] , for hard ellipsoids 
[55] , for the Gay-Berne [56] [57] , and for Yukawa segment 
models [28]. In order to do this one has to map the sys- 
tem at given temperature and density onto the parame- 
ter space. This mapping needs the full direct correlation 
function of the isotropic phase as an input. The latter 
can be obtained either by simulation and liquid-integral 
equations for anisotropic systems [58] or by fundamental- 



measure-theory for anisotropic hard particles [20] . 

In future work the model can be applied and solved for 
more complicated situations such as interfaces between 
two coexisting phases. While the isotropic- nematic inter- 
face has been studied by theory, simulation and experi- 
ment [59] [60] , it would be in particular interesting to get 
structural information about solid-fluid and solid-solid 
interfaces. Then the interfacial tension is anisotropic and 
the degree of anisotropy can be extracted from the phase- 
field-crystal approach. The structure of topological de- 
fects would be another playground where the present ap- 
proach could be applied to [45] [61] . 

Moreover, the extension towards dynamics is straight- 
forward [31] and the numerical implementation can in 
principle be done. We expect a wealth of different dy- 
namical growth effects [62] and novel steady states in 
systems driven by external [63] and internal [30] [64] [65] 
forces. Finally a big challenge is to implement the func- 
tional in three spatial dimensions where there are more 
coupling parameters [32] . 
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